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The recent observation of the dynamical Casimir effect in a modulated superconducting waveguide, 
coronating thirty years of world-wide research, empowered the quantum technology community with 
a powerful tool to create entangled photons on-chip. In this work we show how, going beyond the 
single waveguide paradigm using a scalable array, it is possible to create multipartite nonclassical 
states, with the possibility to control the long-range quantum correlations of the emitted photons. 

In particular, our finite-temperature theory shows how maximally entangled NOON states can be 
engineered in a realistic setup. The results here presented open the way to new kinds of quantum 
fluids of light, arising from modulated vacuum fluctuations in linear systems. 
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I. INTRODUCTION 

Generation of photons out of a perturbed vacuum HHni, usually referred to as dynamical Casimir effect (DCE) 
P^SHTT] is one of the most fascinating predictions of quantum electrodynamics [18] , closely related to the better known 
radiation of Hawking around a black hole [19] . 

The quest for the observation of such an effect has been one of the great scientific endeavors of our times, success¬ 
fully concluded with the unambiguous observation of Casimir radiation emitted out of a modulated superconducting 
microwave resonator isniEi]. Such a groundbreaking result not only definitely vindicated thirty years of study of 
DCE physics, but it opened a new world for quantum technology research, turning DCE from a theoretical quantum 
phenomenon into an useful tool for quantum device engineering [22] |23] . Casimir radiation can be interpreted as a 
parametric amplification of the vacuum, creating pairs of entangled photons out of a perturbed vacuum [24]. This 
makes of DCE an ideal source of entangled photons naturally integrated into superconducting waveguides [254(27] . 

In this work we aim to make a major step forward in this direction, proving how it is possible to harness the existing 
technology of superconducting circuits to both investigate many-body photonic processes and to create nonclassical 
field states with direct applications in quantum communication devices. We propose a scalable lattice architecture of 
modulated coupled waveguides in which it is possible to have a fine control over the long-range photonic correlations 
of the Casimir photons, leading to the possibility to emit on-demand NOON [28] and other nonclassical field states 
and to explore quantum many-body phenomena in coupled linear systems. 


II. RESULTS 


A. The DCE in an array of coupled stripline waveguides 


In this paper we investigate the physics of a one dimensional dynamical Casimir array, that is a one dimensional 
array made of N superconducting open stripline waveguides (SW), each one terminated on one side by a SQUID loop, 
threaded by external flux 4>a(t) [29]. Each terminating SQUID is then coupled to the next one through an additional 
coupling SQUID threaded by an external flux 4>b(t) [30] . 

This arrangement allows us to both modulate the terminating impedance on each SW, and the coupling between 
neighboring SWs. We aim to study how the correlations between Casimir photons coming out of the system through 
the same as well as through different SWs vary as a function of the driving external fluxes. While a priori each SW 
could be driven differently, in this paper we will consider both 4>a(t) and 4>b(t) to be independent from the specific 
SW. As we will see, this arrangement allows us to have a solid control over the behavior of the quantum correlations 
whilst limiting to the minimum the number of independent control knobs. A sketch of the system can be found in 
Fig. 0 It is interesting to notice that arrays of coupled, static quantum circuits (without modulations) have been 
proposed as promising systems in the field of quantum simulations [344(33] . 

Using the standard theory of circuit QED [34l [35], we can identify the dynamical variables of the system as the 
node fluxes at position x along the i-th SW, and their conjugate momenta Px^i (see Eig. |^. The Hamiltonian 
of the system can thus be written as = Hsw + Hsq, where Hsw is the Hamiltonian describing the SWs and Hsq 
the one for the SQUIDs [29] . 
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In Eq. Q, Co and Lq are the characteristic unit length capacitance and inductance of each SW, Ej[4>a(t)] and 
Fj[4>b(t)] are the Josephson energies of the terminating and of the coupling SQUIDs, modulated by the respective 
external fluxes, and Cj is the total capacitance in each SQUID [221 EEl- Notice that, for sake of simplicity, we 
neglected additional smaller capacitive terms in the Hamiltonian due to the presence of the coupling SQUIDs and we 
assumed that the plasma frequency of the SQUIDs far exceeds other characteristic frequencies in the circuit, so that 
oscillations in the phase across each SQUID have small amplitude. 

Our aim is to investigate the correlations among the Casimir photons emitted out of the system through the different 
SWs. We can thus apply an input-output formalism, effectively tracing out all the internal, unobservable degrees of 
freedom. In order to do this we extend the procedure described in Ref. [29], deriving from Hsq a set of boundary 
conditions over the SWs variables. As detailed in the Appendix A, this procedure allows us to obtain a linear algebraic 
system linking bosonic ladder operators for photons getting in and out of the Uth SW, respectively noted and . 
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FIG. 1. Schematic of the array of superconducting waveguides with tunable electric lengths and couplings. The array can have 
open boundary or can be ring shaped. Each SW is terminated by a SQUID and it is coupled to the nearby SW by a second 
SQUID. The SQUIDs impose boundary conditions in the SWs, that can be parametrically tuned by changing the externally 
applied magnetic fluxes. We highlighted in red the thin insulating layer of each Josephson junction. The equivalent circuit 
diagram for the case of two SWs can be found in Fig. 



FIG. 2. Equivalent lumped element circuit diagram for the case of two SWs. The two SWs are characterized by the dynamical 
fluxes and respectively. Each superconducting SW is characterised by its characteristic inductance Lo and capacitance 
Co per unit length. 


Such a system can be easily solved by expressing the coordinates and momenta for each SWs in terms of collective 
coordinates and momenta which diagonalize the Hamiltonian in Eq. 0 m- Input and output operators on the Uth 
SW can then be written as linear functions of the corresponding eigenmode operators as 

N 

= ( 2 ) 

n=l 


where n indexes the eigenmodes and the are real coefficients. 

The SQUID loops provide effective Josephson coupling energies Ej(t) = Ej[^a{t)] and Ej{t) = Fj[^ 5 (t)], tuned 
by the threading magnetic fluxes ^a{t) and with a harmonic time dependence. The resulting time-dependent 
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Josephson energies can thus be expressed as 

Ej{t) =Ao sin((/)) + Mo sin(6>) cos , (3) 

Fj{t) =Ao cos((/)) + Mo cos(6>) cos (ujdt) , 

where, in the perturbative regime we consider, SAq Aq. The angles (j) and 0 describe the ratios between the 
static and variable parts respectively of the terminating and coupling SQUIDs. The resulting output fields are 
pairwise correlated, each one with the mode with angular frequency symmetrical around half the driving frequency 
<^± = ^d/2 ± M. In the degenerate case M = 0 (co’+ = uj- = ojd12)^ the system can finally be diagonalized obtaining 
input-output relations for the normal modes of the system 


Tout _ _Tin _ rr un f 


( 4 ) 


where v is the phase velocity and is a linear function of 5Ao and depends on 9 and 0 (see Eq. § and Appendix 
A). From Eq. we clearly see that, as expected, the time-modulation is responsible for the mixing of creation and 
annihilation operators and thus for the emission of Casimir photons. Notice that the degenerate case, on which we 
will concentrate in the following, is the most robust against thermal noise, since lower energy modes are exponentially 
more affected by thermal noise (See Fig. 6). From Eq. ([^ and Eq. 0 we can finally find the intensity emitted out 
of the Tth SW as 
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Once we have calculated the input-output relations linking input and output channels over the different SWs, we 
are in measure to calculate correlation functions to extract information on the nature of the Casimir photons emitted 
in the circuit. Notice that, while the physically measured quantities are currents and voltages out of each SW, here 
we will present correlation functions involving directly the single mode photon operators calculated from Eq. 0 
and Eq. 0. This allows us to deal with dimensionless quantities and standard quantum optics normal ordered 
correlation functions. Their direct link with correlation functions involving voltage operators is detailed in Appendix 
E. Although intensity detectors (measuring normal order correlation functions) in the microwave frequency range 
are under development [38l |39] , it has been shown that these normal order correlation functions can also be inferred 
by currently used linear detectors [40]. The second order correlation function involving fields from SWs i and j, 
Gpj = ^ vacuum input, can thus be written as 


^5- = (IJ)' E SK SL 

n,m 


(6) 


The aim of the rest of this paper will be to investigate the very rich physics of Eq. 0, and to show how it is possible 
to exploit it to control long-range photonic correlations and to emit nonclassical multiphoton states. 


B. Two coupled stripline waveguides 


We start our analysis from the simplest and analytically solvable case of the DCE array represented in Fig. 
composed of only two coupled SWs, that we will name SWl and SW2 respectively. In this case the transformation 
coefficients in Eq. ([^ are simply cj = = C 2 = l/\/2, = —\l\f2 [37] • We obtain for the intensities. 
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As explained in Appendix B, correlation effects in presence of pairs emission can be better understood considering 
the intra-SW and iter-SW second order normalized correlation functions defined as gf‘^ = such that 

0 < < 1. This normalisation not only allows one to avoid artefacts due to small signals but it also makes the 

second order correlation function independent from the intensity of the amplitudes Aq and 5Aq^ but only dependent 
on (j) and 0. We obtain for the intra-SW and inter-SW second order correlation functions. 
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FIG. 3. (a) Normalized intra-cavity second order correlation function for a system of two coupled SWs as a function of 
9 and 0 for a system of two SWs. (b) Normalized inter-cavity second order correlation function a function of 9 and 0 

for a system of two SWs. (c) Normalised correlation functions g^\ (red continuous line) and g^^\ (blue continuous line) as a 
function of 9 for fixed 0 = 7r/4. These ID plots correspond to the line-cuts indicated by dotted lines in Figs. 3a and b. The 
dashed and dotted lines represent instead correlation functions calculated at T = 25 mK and T = 40 mK respectively, (d) 
Normalized second order correlation functions for a system of five coupled SWs with open boundaries: g^^l (black continuous 
line), g^^l (blue dashed line), and g^^^ (red dotted line). Parameter used in the calculations are Zq = 55 Q, v = 1.2 x 10® m/s, 
Ud — ‘27V X 10.3 X 10^ and, where not otherwise stated, T = 0. Moreover, Aq and SAq were chosen so that the maximum value 
of does not exceed 0.1, at T = 0. 


where 


SLi 


6L2 


1 ^^^0 sin6> 

\27ry I/O Aq sin^ cj) ’ 

/ 00 \ ^ 1 ^^^0 sin 6> + 2 cos 9 
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( 9 ) 


(2) (2) 

and 00 is the flux quantum. When either SLi or 6 L 2 are zero, gH = 9 i 2 ^ that is, after one photon is detected from 
the SWl, there is the same probability to detect at zero delay a second photon from either SWl or SW2. If instead 
6 Ej and 6 Fj are chosen so that 6 L 2 = —SLi, the probability to get two photons from the same SW vanishes. This 
photon blockade effect is usually observed in highly nonlinear systems |4TH44], while the arrays here described work 
in the linear regime. Finally, for 5 L 2 = SLi, the two photons come from the same SW and the probability to find 
one photon in each of the two SWs is zero. Below we will show that this situation corresponds to the realization of a 
NOON stdite [28] 10) = + l^^^)). The system thus works as a switchable source of photons with the unique 

property to dynamically control the correlation properties of the emitted photon pairs. 

The density plots of g^l = g^l and g^l = g^l as a function of 0 and 9 can be found in Figs. 3a and b. We can see 

(2) (2) 

that, for a given 0, changing 9 allows to modulate g^ ( and g[ 2 spectrum of possible values (Figs. 3c and 

d). This is a demonstration of the full quantum control capability offered by the DCE arrays. Note that the two plots 
in Figs. 3a and b are linked by the relation 1 + 2 = 1- Figure 3c shows the line-cuts at fixed 0 = 7r/4 of Figs. 3a 

and b. The angle 0 = 7r/4 corresponds to arrays with a stationary coupling potential equal to the stationary potential 
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of the single unit, hence we find that, in the presence of such a large coupling, situations may exists {0 = arctan (1/4)) 
where photons are bound to their SW and cannot spread in the nearby ones — 0)- 

Since real experiments are carried out at finite temperatures of the order T ^ 20 — 60 mK m, we reported in 
Fig. 3c the inter and intra SW correlation functions at temperatures T = 25 mK and T = 40 mK (details on the 
calculation can be found in Appendix D). These finite temperature calculations demonstrate that the predicted effects 
are solid and can thus be observed in real experiments under the same conditions used to demonstrate the DCE. 

The obtained intensity correlation functions exhibit interesting quantum features. In a two-mode system quantum 
correlations may exist which violate classical inequalities. If the two modes are symmetric, as in the present case, 
according to the Cauchy-Schwarz inequality, two fluctuating classical fields (described by a positive Glauber-Sudarshan 
P function) satisfy the following inequality: 


VI,2 — VI,1 


( 10 ) 


For (j) values such that arctan(—2) ^ O.GStt < 0 < tt, this inequality is violated (see Fig. 2c) hence the system displays 
nonclassical correlations. At angles where 6 Li — 6 L 2 = 0 e.g., 0 = 7 r /4 and 0 = arctan(—1/5), there is the maximum 
violation: no photon pairs in a single SW = g ^2 — 0) perfect inter-SW pair correlation g ^2 — 1- 


C. Beyond two coupled waveguides 


Figure 3d describes normalized second order correlation functions obtained for an open DCE array composed by 
5 SWs. The panel has been obtained for 0 = 7 r /4 and displays ^ 3^3 (black continuous line), g^\ (blue dashed), and 
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FIG. 4. Normalized second order coherence function for ring chain with 31 SWs and ^ — 'k (a) 3D representation of 

in function of site j = —5, ..,+5 and parameter Q. The blue line-cuts focus on three different scenario for (b) Q — 2.78, 
(c) 6> = 2.56, (d) 6> = 0.17. 


g^^\ (red dotted) as a function of d. We note that the intra-SW correlation remain almost the same as that for the 

( 2 ) ( 2 ) 

2 SWs array, and there are regions where both the inter-SW correlations g^ 2 g^ ( are larger than the intra-SW 

one ^ 33 . It is surprising to see that there are 0 values where it is more probable to find pairs in more distant SWs 
rather than in the same SW or in adjacent SWs. The generated photon pairs show indeed an effective photon-photon 
interaction that can be tuned from attractive to repulsive simply by acting on the ratio between the modulation 
amplitudes. However, as discussed above, the possibility to drastically change the inter- and intra-SWs quantum 
correlations is actually due to quantum interference effects. We have also investigated DGE arrays with a larger 
number of units. Figure 4 shows results for a ring chain with 31 SWs. In Fig. 4 a we plot the correlation functions 
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for 0 = 7 r /4 as a function of 0 and j. As also shown by the line-cuts obtained for three specific angles 6 >, the 
spatial photon correlations can be controlled to present very different scenarios. Assuming that a photon escapes from 
the i-th SW, the probability of a second photon from any different SW vanishes for 0 = 0.17 (Fig. 4d). This scenario 
corresponds approximatively (because g\ ■ is sligthly lower than 1 ) to the highly entangled multimode AO 0 A state 
I'lp) = :^(|2, 0,..., 0) + |0,..., 2,..., 0) + • • • + |0,..., 0, 2 )). On the contrary, for 0 = 2.56 (Fig. 4c), the probability 
to find both photons in the same SW is zero. Fig. 4 a shows a still different situation where, for a photon in the j 
unit, the second one is delocalized around the nearest 8 units. 


D. Entanglement 


The above results provide clear indications that, at specific highly entangled two-photon multimode NOON 
states involving all the array units can be realized. In order to confirm this prediction, we have calculated the two- 
photon array density matrix for two coupled SWs. The matrix elements can be calculated in terms of expectation 
values involving SWs output creation and annihilation operators (see Appendix C). Since all the calculations here 
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FIG. 5. The von Neumann entropy Ej\f (black continuos line) for the reduced density matrix of a system of two coupled 
waveguides as function of 0 and the fidelity JToon of the AOOA state calculated at T = 0 K (blue dashed line), T = 50 mK 
(red dotted line), and T = 60 mK (green dashed dotted line). The displayed results have been calculated using 0 = 7r/4. 


presented have been obtained perturbatively, in the limit of small modulation amplitudes, in the derivation of the 
density matrix we considered only those events where at least one photon was generated (post-selection). In this case, 
each subsystems is constituted by three states corresponding to 0, 1, 2 photons. In the zero temperature limit the 
resulting density matrices describe pure states. To quantify the entanglement we calculate the von Neumann entropy 
logsP/c of iFe reduced density matrix with one of the two SWs traced out, where pk are the eigenvalues 
of the reduced density matrix. We considered a 3-basis logarithm so that the maximally entangled state corresponds 
to E^ = 1 . While the entropy Ej^^ of the total system is null at T = 0, indicating that the state is pure, we can see 
from Fig. 5 that the entropy of the reduced density matrix is greater than zero for all values of 0 except 6 > = 2.94. It 
means that {ip) is an entangled state. As shown in Fig. 5, for 0 = 0.038 and 0 = 1.305, the entanglement is maximum. 
This situation corresponds to the generation of the maximally entangled state, 

|V’) = T(|i,i) + |2,0) + |0,2)) , (11) 


where the first (second) entry in the ket indicates the photon number in the first (second) SW. When 0 approaches 
the condition SLi = —SL 2 ^ so that no photon pairs can be generated in the same SW, the entanglement reduces to 
zero. However, as certified by the violation of the Cauchy-Schwarz inequality in Eq. (10), also in this case the system 
exhibits highly nonclassical features. 

Finally we calculated the fidelity of the AOOA state Thoon = ^/{Pnoon), where pnoon = |'0noon)('0noon| with |'0noon) = 
^(|2,0) + |0,2 )). It results (see Fig. 5) that for 0 = arctan(—5) corresponding to the condition 6 Li = (iI/ 2 , the 
fidelity becomes unitary, demonstrating that a two photon NOON state is actually obtained. Figure 5 also reports 
finite temperature calculations at T = 50 mK and T = 60 mK showing that entanglement results are solid against 
temperature and the generated AOOA states can be experimentally observed. 





III. CONCLUSION 


In this paper we studied the physics of an array of SWs in which both the terminating impedances and the inter- 
SW couplings can be independently modulated in time. We discovered that in this way it is possible to control the 
long-range quantum correlations of the emitted photons. In particular, we found it is possible to swipe the second 
order correlation functions between photons, in the same or in different SWs, over their full range, by tailoring the 
modulations. We also proved how entangled multipartite NOON states can be efficiently emitted by the system. 

The results here presented open the way to new kinds of quantum fluids of light [45], where effective photon- 
photon interactions arise and can be actively controlled in linear systems. In these systems photons needs not to be 
externally pumped [32l|33l|46], but they are generated through quantum vacuum stimulation, so that their correlations 
can be nonclassical to the highest degree. Moreover DCE arrays, scalable and versatile by design, naturally allow for 
control, manipulations via external magnetic fields, and measurements of individual lattice sites. For example, the 
open waveguides can be replaced by closed resonators and defects can be easily introduced simply by changing the 
modulation amplitude of one or more coupling or terminating SQUIDs. The present work is only a first investigation 
of the potentialities of DCE arrays. Since the DCE turns vacuum fluctuations into real, observable particles, DCE 
arrays constitute a test bed for the study of quantum-vacuum correlations in many-body quantum systems. Future 
works will have to further explore their possibilities by exploiting their intrinsic design flexibility to produce novel 
kind of nonclassical many-photon states with long-range correlations when the weak-modulation hypothesis fails. 
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Appendix A: Output field operators for DCE arrays 


In the continuum limit {Ax 0) [29|, the flux coordinate of the i-th SW can be expressed as a continuous function 
of the position: 

In order to derive the output fields, it is useful to find the normal modes of the system transforming coordinates 
and momenta according to and Pi{x) = c^P^(x), such that Hs can be expressed in terms of 

independent oscillators. 


Hs 


{t) $0,n + 2 E 


Cj ’ 


(Al) 


where the Kn{t) are linear combination of Ej{t) and Fj{t). Equation (Al) shows that the whole system can be 
regarded as a collection of independent parametric oscillators. By writing the Heisenberg equations of motion under 
the influence of the total Hamiltonian P, and taking the continuum limit, the following boundary condition (at x = 0) 
can be obtained [29|, 


Cj^n (0, t) + An (t) (0, t) 


1 d^r. 


I/O dx 


= 0 . 


(A2) 


a :;=0 


The SW coordinates as well as the eigenmodes collective coordinates can be expanded in terms of input 

and output creation and annihilation operators [2l]. For example the collective coordinates can be expanded as 


fm r 
V dTT I 




bj," (w) + b^* (w) + h.c. 


(A3) 


where Zq = ^JLQ|Co is the characteristic impedance and = co/v, being v the SW phase velocity. The multimode 

input and output photon operators obey bosonic commutation relations = 5n,m 

We now consider applied magnetic fluxes such that Ej{t) = Aq sin(0) sm{0) cos f^t) and Ej( t) = Aq cos(0) + 
^Ao cos(0) cos {ujdt)^ with weak modulation amplitudes at frequency ujd- Inserting Eq. ( |A3[ ) into Eq. ( |A2[ ) and solving 
perturbatively we can derive the normal mode output fields in terms of the input normal mode operators [29] . 

t>r*(w) = + Sn{uJ,UJd - Uj)h'^'‘{uJd - w) , 


(A4) 
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where 


SL 


Sn {uj' ,uj") = —I — -Vu/u/'Q {oo') 0 {oo "). 


(A5) 


In Eq. (A5) 0(co’) is the Heaviside function and 6 Ln = ((/)o/27r)^ SAn/ (I/qAq^), where we have defined A^ = Aon + ^^A^, 
with Aon = hm^Ao^o An- 

In order to derive adimensional photon-photon correlation functions it is useful to define single frequency output and 
input photon operators corresponding to the creation or destruction of photons within a small frequency bandwidth 
A around a central frequency ujj [24] 


bniuij) = {1/VA) / duhn(ui), (A6) 

obeying standard bosonic commutation relations (uif)] = 6n,mS: In the main body of the paper we 

concentrate our attention to the degenerate case with the Casimir photon pairs both at frequency uJd/2 and thus, for 
the sake of simplicity, we define bn{oJd/‘^) = ^n- 


Appendix B: Normalized correlation functions 

The second order correlation function for a generic field a is usually normalized over the squared population: 

’ leading to = 1 for a perfectly coherent field. This choice of the normalization is 
nevertheless not well suited to the case of weak processes creating pairs of excitations. The problem can be easily 
understood considering the state in which parametric process creates a pair of excitations with low probability a‘^ 1: 

\^) = |0) +a|2), where |n) is the state with n quanta (such a state is normalized to Direct calculation shows 

that = (2(a)“^, that is, the correlation diverges for a vanishing parametric process. A more physically meaningful 
way to normalize the second order correlation for these processes, that is the one we employ in the present paper, is 
g^‘^^ = ; f^at is bounded to unity for states with up to two photons. 


Appendix C: Derivation of the density matrix 

In order to calculate the von Neumann entropy and the fidelity of NOON states to quantify the entanglement, 
we derived the density matrix p for the entire system. For a weak driving, a perturbative calculation describing 
the generation of a single photon pair is adeguate, and the Hilbert space for the two SWs system consists of the 
tensor product of two qutrits (three-level systems). Each SW can have only 0, 1, 2 photons. The matrix elements 
Pnm^n'm' = {\n'm'){nm\) (with n^m^n'^m' = 0,1,2), can be expressed in terms of one and two photon correlation 
functions already used in the text. For example (|20)(11|) = (l/v^)((a|)^aia 2 ). Since all the calculations here 
presented have been obtained perturbatively in the limit of small modulation amplitudes, the matrix element poo,oo 
is the larger one and we get rid of it by post-selecting only those events where photons are emitted. 


Appendix D: Finite temperature 


Finite temperature calculations, displayed in Fig. 3c and 5, can be carried out as the zero-temperature ones, with 
the only difference that the normal order mean values of the input photon numbers (at frequency ujd/^^) acquire the 
thermal equilibrium values where Nt = -g thermal photon occupation 

of the input-field mode with frequency T is the temperature, and ks the Boltzmann constant. Higher order 

input correlation functions can be simply derived employing ordinary Gaussian factorization and using the thermal 
equilibrium result {b'^n^m) — ^ Hi]- We obtain for the first order correlation function. 


GP=y:{<) 


Nr + il + Nr) 


(Dl) 
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The second-order correlation function is analogously given by 


= E E <<cL4nC^ Nt + If SK5L, 

n m 


(c^ + C* C* 1 

n m 


Nr+y^y 5Ll {Nt + 1) 


Nr+y^y 5Ll {Nt + 1) 


(D2) 


Appendix E: Voltage-based correlation functions 

In microwave circuits, the physically measured fields are the current and voltage along the line. The voltage operator 
for each SW can be obtained directly from the corresponding flux field through the relation (x, t). The 

signal recorded by a photon intensity detector in the i-th SW corresponds to the normal order correlation function 
Gf^ (r) = {Vr{T)V+{Q)), where y+ and W are the positive and negative frequency components of the i-th SW 
output voltage + Vy ~, which can be expressed in terms of annihilation and creation operators respectively. 

Although intensity detectors (measuring normal order correlation functions) in the microwave frequency range are 
under development [38], it has been shown that these normal order correlation functions can also be inferred by 
currently used linear detectors m- Following this procedure, higher order correlations involving voltages also from 
different SWs can be obtained. In the main text, in order to present adimensional normalised correlation functions, we 
used single mode correlation functions involving photon operators for a small bandwidth around a specific frequency 
Here we present voltage-based correlation functions involving fields with a broader spectral range. By using 
the output photon operators derived in the Methods section, it is possible to derive the first order correlation function 


+ (r) = ^ 

where Zq = y/LoJCo is the characteristic impedance. We obtain 


POO POO 

Cffr) = ^ / / (w') a?"* {co")}e<^'-^'> , 

Jo Jo 


(El) 


= r . (E2) 

Jo 

The photon flux spectral density in the output field, defined by [29] 

POO 

nrfco) = / da;'(art^)ar(^')) (E3) 

Jo 

can thus be expressed as 

1^ • (E4) 

n 

Figure 6a displays the photon flux spectral density as a function of the normalized detection frequency {ujd is the 
modulation frequency) calculated for a system of two SWs for T = 0 (black solid line), T = 25 mK (blue solid line) 
and in the absence of modulation at T = 25 mK (red dashed line). Parameters are provided in the figure caption. 

The second order correlation function, defined as = (K-“(r)i/“(r)l/^(0)V-^(0)) can analogously be written 

as 

G'!f(r)=(^) ^44 c4c)„J„(t)C(t), (E5) 

^ ' n,m 

where 

r^d 

^n{r) = / dcjy/cj {uJd - ^)Sn (^, “ ^) . (E6) 

Jo 

Figure 6 displays G^^^(r)/G^^^(0) and G^2^(r)/G^^^(0) as a function of the time delay for a system of two SWs. The 
figure shows a decay of these correlations in the absence of rapid oscillations and without significant modifications of 
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FIG. 6. The parameters used in the calculations are Zo = 55 u = 1.2 x 10® m/s, = 27r x 10.3 x 10^. Aq and SAq were 
chosen so that the maximum value of does not exceed 0.1 at T = 0. (a) Photon flux spectral density as 

a function of the normalized detection frequency calculated for a system of two SWs for T = 0 (black dotted line), T = 25 
mK (blue solid line) and in the absence of modulation at T = 25 mK (red dashed line), (b) Second order correlation function 
(black line) and (blue solid line) as a function of the time delay for a system of two SWs. 



e 


FIG. 7. Zero-delay normalized correlation functions (A)/y and Gl' 2 ^( 0 )/^GyGy as a function of the parameter 

0 are plotted for a system of two SWs at zero temperature. 


the ratio G^^\r)/G^^\r) when time increases, analogous to what is observed in a single waveguide [29]. These results 
show that the zero delay correlation functions presented in the main text are observable in real experiments where 
nonzero time-windows are required. 

The zero-delay normalized correlation functions {G )/and (0)/as a function of the 
parameter 0 are reported in Fig. 7. The figure shows that a behavior analogous to the single frequency case showed in 
the main text (see Fig. 3c) can also be observed in broad-band voltage measurement. The only significant difference 
is that these normalized correlation functions are no more bounded by 1. 
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